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SUMMARY 


In  this  report  we  detail  the  development  of  a  new  regional  3-D  tomographic  P-wave 
velocity  model  (WINPAK3D)  of  a  region  in  Southern  Asia  centered  on  Pakistan.  Our 
primary  goal  in  developing  a  3-D  model  of  the  crust  and  upper  mantle  in  this  region  is  to 
improve  regional  seismic  event  location.  The  initial  model  for  the  region  was  developed 
by  integrating  the  results  of  more  than  60  previous  studies  related  to  crustal  and  upper 
mantle  velocity  structure.  We  refined  the  model  by  applying  a  joint  velocity  tomography 
and  event  location  procedure  to  a  dataset  of  earthquakes  in  the  region.  The  model  was 
iteratively  updated  using  a  nonlinear,  conjugate-gradients  technique  that  adjusts  the 
velocity  model  to  minimize  the  misfit  between  calculated  and  observed  travel  times  from 
multiple  stations  and  events,  subject  to  smoothness  constraints.  We  compute  the  travel 
times  and  their  sensitivities  to  the  velocity  structure  with  an  extension  of  the  3-D  Podvin- 
Lecomte  (1991)  method.  Because  the  event  locations  are  not  fixed,  we  also  relocate  all 
events  using  a  3-D  grid  search  location  method  after  each  update  of  the  3-D  velocity 
model.  The  data  set  consisted  of  a  suite  of 535  events  containing  over  6,600  arrivals 
from  the  Engdahl  et  al  (1998)  database.  We  performed  several  iterations  of  the 
nonlinear  algorithm  to  invert  for  the  Pn  velocity  as  a  function  of  latitude  and  longitude 
and  then  imposed  a  fixed  rule  for  extrapolating  this  velocity  into  the  upper  mantle. 

Future  extensions  of  the  method  will  allow  for  more  flexible  changes  in  the  upper  crust 
and  mantle  portions  of  the  model.  Our  results  show  that  the  new  model  better  fits  the 
data  compared  to  both  the  initial  model  and  the  global  1-D IASP91  model.  The  root 
mean  square  (RMS)  error  for  the  the  updated  3-D  model  is  1.81,  compared  to  2.01  for  the 
initial  3-D  model  ami  2.42  for  the  IASP9I  model. 


INTRODUCTION 

An  accurate  estimate  of  the  location  of  a  seismic  event  is  particularly  important  for 
monitoring  potential  nuclear  explosions.  Determining  hypocenters  for  small  seismic 
events  (mb  <  4.0)  with  high  accuracy  is  difficult  because  of  limited  recordings  and 
complicated  crustal  structure  at  regional  distances  that  is  not  accounted  for  in  the 
standard  global  1-D  models  such  as  the  IASP91  model  (Kennett  and  Engdahl,  1991). 
Precise  regional  location  of  seismic  events  requires  a  velocity  model  that  accurately 
represents  the  real  Earth  structure,  as  systematic  biases  caused  by  unmodeled  Earth 
structure  are  known  to  play  an  important  role  in  earthquake  location  errors  (e.g.,  Douglas, 
1967;  Dewey,  1972;  Engdahl  and  Lee,  1976;  Jordan  and  Sverdtup,  1981;  Pavlis,  1992). 
Regional  3-D  models  that  better  represent  the  true  Earth  structure  can  be  used  to  compute 
accurate  travel  times  of  regional  seismic  phases  such  as  Pn,  Pg,  Sn,  and  Lg.  Travel  times 
calculated  using  3-D  models  can  that  be  used  to  develop  source  specific  station 
corrections  (SSSC’s)  that  can  be  implemented  by  monitoring  organizations  to  provide 
improved  regional  event  locations.  Our  goal  in  this  research  project  was  to  develop  a 
regional  model  of  the  crust  and  upper  mantle  for  the  India-Pakistan  region  that  will 
improve  event  location  accuracy.  We  built  our  model  for  the  India-Pakistan  region  by 
updating  a  detailed  preliminary  model  through  a  joint,  nonlinear  velocity  tomography  and 
hypocenter  relocation  technique.  In  the  following  sections  we  describe  first  the  inversion 
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algorithm  followed  by  the  application  of  the  algorithm  to  a  set  of  earthquake  data 
surrounding  a  region  centered  on  Pakistan  contained  inside  25-40°  N  and  60-75°  E. 

JOINT  VELOCITY  TOMOGRAPHY/EVENT  LOCATION 
ALGORITHM  DEVELOPMENT 

The  development  of  high-resolution,  3-D  velocity  models  using  regional  earthquake  data 
is  complicated  by  the  dependence  of  the  model  on  mislocations  of  the  earthquake 
hypocenters.  Travel  time  measurements  from  earthquakes  depend  both  on  the  earthquake 
locations  and  the  earth's  velocity  structure.  Hypocenter  estimates  are  biased  by  errors  in 
the  velocity  model.  Conversely,  tomographic  images  of  the  velocity  structure  derived 
from  earthquake  arrival  times  are  affected  by  errors  in  the  earthquake  hypocenters.  It  has 
long  been  recognized  that  earthquake  location  and  velocity  imaging  are  coupled  inverse 
problems.  The  problems  separate  only  when  sufficient  prior  data  exists  to  constrain  one 
set  of  parameters  or  the  other,  or  in  special  geometries.  Our  study  in  Pakistan  was 
motivated  by  the  lack  of  prior  information  needed  to  decouple  the  location  and  velocity 
imaging  problems.  Furthermore,  we  wanted  to  use  the  wealth  of  information  contained 
in  local  and  regional  earthquake  data,  whose  source-receiver  geometry  prevents  the 
separation  of  the  problems.  Our  inversion  algorithm  jointly  solves  the  nonlinear  problem 
through  the  application  of  a  conjugate  gradient  technique,  which  iterates  over  linear 
inversion  steps  that  include  updates  of  the  hypocenters,  velocities  and  ray  paths.  The 
technique  explicitly  addresses  the  coupling  between  the  hypocenters  and  the  velocity 
structure  by  computationally  breaking  down  the  large  matrix  that  must  be  inverted  for 
velocities  and  locations.  This  computational  technique  results  in  two  separate  smaller 
matrices  which  may  be  inverted  separately,  but  still  solve  the  simultaneous  problem 
(Spencer  and  Gubbins,  1980;  Rodi  etal.,  1981). 

Figure  1  is  a  ample  flowchart  of  the  algorithm  we  used  to  develop  our  new  3-D  velocity 
model  of  the  India-Pakistan  region.  There  are  three  major  components  involved  in  our 
joint  tomography/location  procedure:  1)  3-D  ray  tracing  to  predict  first  arrival  times  using 
an  enhanced  version  of  the  Podvin-Lecomte  (P-L)  method  (1991);  2)  a  3-D  grid  search 
location  algorithm  (GSEL)  from  Rodi  and  ToksOz  (2000)  to  relocate  events  inside  the 
appropriate  velocity  model;  and  3)  a  linear  conjugate  gradient  inversion  algorithm  to 
produce  the  updated  velocity  model  inside  each  iteration  of  the  overall  nonlinear 
algorithm.  In  the  next  few  subsections  we  describe  the  velocity  model  parameterization 
and  each  of  the  three  major  algorithmic  components. 
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Figure  1:  Simplified flowchart  of  the  nonlinear  joint  tomography  and  location  procedure 
used  to  develop  the  3-D  model  of  the  region  in  and  around  Pakistan. 


Velocity  Model  Parameterization 

We  parameterize  our  3-D  model  of  the  Pakistan  region  in  terms  of  a  velocity  vs.  depth 
profile  at  each  point  on  a  geographic  grid,  which  is  sampled  uniformly  in  latitude  and 
longitude.  At  each  geographic  grid-point,  the  velocity  profile  is  given  as  velocity/depth 
pairs  at  nodes  ranging  from  sea  level  to  a  depth  of 760  km.  Discontinuities  in  velocity  are 
allowed  at  the  ocean  bottom,  Moho  and  the  major  mantle  discontinuities  (410  and  660 
km  depth  in  the  IASP91  model).  Using  this  parameterization  also  reduces  the  number  of 
model  parameters  to  be  solved  for  in  the  inversion.  One  complication  in  using  this 
geographically-defined  velocity  model  is  that  our  current  implementation  of  the  Podvin- 
Lecomte  algorithm  solves  the  eikonai  equation  in  Cartesian  coordinates  for  a  flat-earth 
model.  However,  we  have  developed  the  algorithms  for  accurately  mapping  our 
geographic  velocity  model  to  a  Cartesian  block  model  required  by  the  P-L  algorithm  and 
for  mapping  the  3-D  Cartesian  travel  time  grids  and  related  ray  path  sensitivities  back  to 
geographic  grids  for  use  in  the  inversion. 


Inversion  Method 


We  formulate  our  inversion  approach  as  follows:  the  unknowns  are  a  vector  m  containing 
the  velocity  model  parameters  to  be  estimated  and  the  hypocenters  (x/)  and  origin  times 
(tj)  of  M  events:  (%,  tj),j  -  1,...,  M.  The  data  are  arrival  times,  d,h  from  each  event  to  a 
subset  of  N  stations  indexed  as  i  =  1,...,  N.  The  data  and  unknowns  are  related  by 

d*='/+:T/(x,;m)+e#,  0) 

where  ei}-  is  the  error  in  dy  and  7}  is  a  function  that  predicts  the  travel  time  to  station  /  from 
an  event  hypocenter  x,.  This  function  depends  on  the  model  parameter  vector  m.  Our 
joint  inversion  criterion  is  to  minimize  an  objective  fimction  of  the  form 

¥(m, x, , /„ . xM , tM )  =  £ |  dv  -  tJ  +  T,  (\j ;  m)  |2  !a\ +  r  \  Lm  |2  (2) 

9 

with  respect  to  all  the  unknowns.  Here,  O/y  is  the  standard  deviation  of  ey.  The  second 
term  of  \| /  imposes  a  smoothness  constraint  on  the  velocity  model,  where  L  is  a 
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regularization  operator  and  r  a  regularization  parameter  The  operator  L  is  a  differencing 
operator  that  minimizes  the  spatial  derivatives  of  the  model  velocity.  The  parameter  t 
determines  the  degree  of  model  smoothness. 

We  perform  the  minimization  of' ¥  numerically  using  a  combination  of  nonlinear 
conjugate  gradient  (NLCG)  and  grid-search  techniques.  The  NLCG  technique  is  used  to 
update  the  model  m  iteratively  along  a  sequence  of  computed  search  directions.  Inside 
each  NLCG  iteration,  three  processes  occur:  3-D  ray  tracing  to  predict  travel  times,  a  grid 
search  location  method  to  minimize  \p  with  respect  to  all  the  event  hypocenters  (x,)  and 
origin  times  (//)  with  the  model  fixed,  and  a  linear  conjugate  gradients  method  to  update 
the  velocity  model  m.  The  grid  search  for  a  given  event  is  performed  within  a  specified 
epicentral  radius  and  depth  range  from  its  initial  location,  allowing  us  to  handle  events  of 
varying  ground-truth  levels  (e.g.,  GTO,  GT5,  GT15).  The  linear  conjugate  gradient 
method  iteratively  computes  an  optimal  update  to  the  current  3-D  velocity  model  using 
the  appropriate  travel  time  tables  and  hypocenters  for  that  model.  This  optimized  update 
is  then  used  as  a  search  direction  in  the  NLCG  inversion. 

We  note  that  our  joint  inversion  algorithm  is  fully  nonlinear  with  respect  to  both  the 
velocity  model  and  event  locations  since  travel  time  modeling  and  event  relocation  are 
performed  for  each  update  of  m.  We  also  note  that  we  are  currently  solving  only  for  the 
Pn  velocity  as  a  function  of  latitude  and  longitude  and  imposing  a  fixed  rule  for 
extrapolating  this  velocity  into  the  upper  mantle.  Future  extensions  of  our  approach  will 
implement  more  general  model  parameterizations  to  allow  for  changes  in  Moho  depth 
and  more  general  variations  in  upper  mantle  velocity. 


Travel  Time  Calculation  Using  Podvin-Lecomte  Rav  Tracing 

To  determine  the  values  of  the  travel  time  functions  Tt,  we  evaluate  Tt  (x;  m)  for  a  fixed 
hypocenter  by  interpolating  a  travel  time  table  stored  on  a  3-D  hypocenter  grid.  This  grid 
is  created  by  applying  the  Podvin-Lecomte  (P-L)  finite-difference  travel  time  algorithm 
(Lomax,  1999;  Podvin  and  Lecomte,  1991)  to  the  earth  model  defined  by  m,  using  the 
location  of  the  rth  station  as  the  “source”  in  the  calculation.  The  Podvin-Lecomte  method 
solves  the  eikonal  equation  in  a  3-D  medium  using  a  finite-difference  approximation.  It 
can  accurately  model  different  propagation  modes,  such  as  transmitted  and  diffracted 
body  waves  or  head  waves.  It  estimates  accurate  travel  times  in  the  presence  of  severe, 
arbitrarily-shaped  velocity  contrasts,  as  occur  across  the  Moho  discontinuity.  This  is  an 
improvement  over  other  similar  methods  (Vidale,  1988, 1990;  Moser,  1991),  which  can 
encounter  serious  difficulties  in  the  presence  of  sharp  first-order  contrasts.  The  model  is 
discretized  on  an  equally  spaced  grid  comprised  of  constant  velocity  cells.  Multiple 
arrivals  (transmitted,  diffracted,  and  head  waves)  are  calculated  at  each  grid  node  and  the 
first  arrival  time  is  chosen.  The  time  t  at  the  current  node  is  a  function  of  the  times  t„  at 
some  (3  or  fewer)  of  the  neighboring  nodes  and  the  slowness,  s,  in  the  cell  traversed  by 
the  wavefront  to  reach  this  node.  That  is,  t=t  (t„,  s).  This  method  of  travel  time 
computation  produces  a  full  grid  of  travel  times  considerably  faster  than  two-point  ray 
tracing,  and  the  sources  and  receivers  can  be  located  anywhere  within  the  model.  The 
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Podvin-Lecomte  computations  are  output  in  the  form  of  3-D  travel  time  tables,  one  for 
each  station  and  seismic  phase,  which  can  be  that  used  by  a  grid-search  event  location 
algorithm  in  lieu  of  global  travel  tables  such  as  IASP91 . 

We  have  extended  the  P-L  algorithm  to  compute  the  sensitivities  of  travel  times  to  cell 
velocities,  as  required  for  the  inversion.  To  calculate  the  ray  path,  we  save  the  node 
pattern  (“stencil”)  at  each  step  through  the  model  as  we  perform  the  normal  P-L  forward 
travel  time  calculation  to  predict  arrival  times.  This  stencil  indicates  which  of  the 
neighboring  nodes  were  used  to  calculate  the  minimum  time  at  the  current  node.  The 
stencils  can  be  used  to  reconstruct  a  path  from  any  node  of  the  grid  to  the  source.  The  ray 
tracing  is  accomplished  by  identifying  all  of  the  grid  nodes  and  the  cells  (slownesses)  that 
contribute  to  the  calculation  of  the  time  at  the  receiver.  As  the  wave  front  propagates 
away  from  the  source,  more  nodes  (and  cells)  are  involved  in  the  travel  time  calculation 
at  each  step.  After  the  midpoint  of  the  ray  path,  the  propagation  region  narrows  until  it 
readies  the  one  node  at  the  source  location.  The  sensitivity  of  the  travel  time  to  the 
slowness,  dt/ds,  is  calculated  at  each  grid  node  of  the  "ray”  for  the  last  cell  traversed  by 
the  wavefront  to  reach  that  node.  The  weight  of  each  neighboring  node  in  the  calculation 
of  the  time  at  the  current  node  ( dt/dt„ )  determines  the  weight  of  the  subsequent  node-to- 
source  subpath  in  the  total  travel  time  calculation  for  the  ray.  The  sensitivities  along  each 
subpath  are  then  weighted  by  this  term. 

Figure  2  shows  XZ,  YZ  and  XY  projections  of  the  ray  sensitivity  matrix  at  station  ASH 
for  an  event  on  the  Afghanistan-Pakistan  border.  Note  that  the  ray  paths  are  not  straight 
lines  but  are  instead  "tubes"  of  the  region  that  contribute  to  the  calculation  of  the  shortest 
travel  time  between  two  points.  Figure  3  shows  the  sensitivities  of  the  travel  times 
calculated  by  P-L  to  cell  slownesses  for  paths  originating  at  station  DSH  in  Tajikistan. 

The  rays  exhibit  strong  sensitivity  to  the  model  near  the  Moho  boundary  (presumably  in 
the  Pn  velocity  depth  range).  This  suggests  that  updates  to  the  Prt  velocities  as  well  as 
the  Moho  depth  will  provide  the  most  significant  improvement  to  the  velocity  model. 

After  the  sensitivities  are  calculated  for  all  of  the  station  Cartesian  grids,  they  are  then 
mapped  to  a  geographic  grid.  Following  this  procedure  Pn  sensitivities  are  extracted  to 
form  the  kernel  matrix  for  the  tomographic  inversion. 

The  ray  tracing  algorithm  and  sensitivity  calculation  were  tested  for  precision  by 
comparing  the  travel  time  computed  as  the  sum  of  the  weighted  sensitivity-slowness 
products  to  the  forward  P-L  calculated  times.  For  rays  with  105  nodes,  the  difference 
between  travel  times  calculated  by  these  two  methods  is  on  the  order  of  10'2  s  when  the 
calculation  is  done  in  single  precision.  Error  of  this  magnitude  for  the  number  of  nodes 
in  the  ray  can  be  accounted  for  by  the  level  of  precision  in  the  calculation,  demonstrating 
that  the  ray  tracing  technique  is  accurately  tracing  the  minimum  time  P-L  ray  path. 
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Figure  2:  XZ  (top),  YZ  (middle)  andXY  (bottom)  projections  of  the  ray  path  sensitivity 
matrix.  This  ray  corresponds  to  an  event  near  the  Afghanistan/Pakistan  border  recorded 
at  station  ASH. 
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Figure  3:  Ray  path  sensitivities  calculated  using  the  extended  P-L  algorithm,  (a)  Station 
DSH  and  events  for  which  sensitivities  were  calculated  shown  with  straight  line 
approximations  of  the  ray  paths,  (b),  (c),  and  (d)  Sensitivities  from  50  km,  60  km  and  70 
km  depth  slices  calculated  for  station  DSH  and  events  shown  in  panel  (a).  Note  that  the 
sensitivities  are  unitless.  The  sensitivities  for  all  rays  that  encounter  a  cell  are  summed, 
producing  scales  that  can  range  from  zero  upward  Deeper  slices  through  the  sensitivity 
matrix  naturally  have  smaller  scale  ranges,  since  the  rays  spread  out  as  they  propagate 
away  from  the  station. 


3-D  Grid  Search  Event  Location 

Our  approach  to  seismic  event  location  is  based  on  a  maximum-likelihood  formulation. 
We  define  an  optimal  location  for  a  seismic  event  to  be  that  which  maximizes  a 
likelihood  function,  constructed  on  the  basis  of  an  assumed  statistical  model  of  errors  in 
the  seismic  data.  Confidence  regions  are  defined  in  terms  of  hypothesis  tests  using 
likelihood  ratios  as  the  test  statistics. 

The  origin  parameters  of  a  seismic  event  are  a  3-D  hypocenter  vector  x  and  an  origin  time 
t.  Given  n  arrival  time  data,  di,  the  location  problem  for  the  event  may  be  expressed  as 

dt  =  t  +  Tfx) +e„i  =  (3) 


11 


where  Tt  is  a  travel  time  function  (travel  time  table)  for  the  rth  datum,  and  e,  is  the 
observational  (picking)  error.  The  index  /  counts  over  both  stations  and  phase  types  (P, 
S,  etc.). 

Our  grid  search  event  location  (GSEL)  algorithm  assumes  the  picking  errors  are 
statistically  independent  and  have  an  exponential  power  distribution,  whose  probability 
density  function  (p.d.f.)  is  given  by 


/feK]  = 
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lie 

> 

P 

► 

K(P)*,  “P 

p\o 

-  ► 
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where  p  >  1  and  K(p)  -  2p1/p  r(J +I/p).  A  Gaussian  p.d.f.  results  with  p= 2,  and  an 
exponential  p.d.f.  with p=\.  We  assume  the  standard  errors,  <r„  are  known  within  a 
constant  factor  and  write 


G;  =  avi  (5) 

where  the  vt  are  known  but  a  is  not. 


The  unknowns  in  the  single-event  location  problem  are  the  origin  parameters,  x  and  t,  and 
the  standard  error  o.  We  allow  a  priori  constraints  on  o,  origin  depth  z  and  origin  time  t 
in  the  form  of  upper  and  lower  bounds.  The  event  epicenter  can  also  be  restricted  to  be 
within  a  specified  distance  of  a  given  geographic  location. 

Considered  as  a  function  of  the  unknown  parameters,  the  joint  p.d.f.  of  the  data  defines  a 
likelihood  function,  which  we  denote  Z.(x,f,o  |  d).  Our  assumptions  imply  this  function  is 
given  by 

-  log  £(x, /,  a  |  d)  =  Y  logo,  +n\ogK(p)+n\o%G  + — ^- *?(*,/ |d)  (6) 

m  P<r p 

where  *F  is  a  data  misfit  function  given  by 

¥(m  |d ) = JtoY  (?) 

i=l  / 

The  optimal  estimates  of  the  unknown  parameters  are  those  maximizing  L,  subject  to  the 
prior  constraints.  For  x  and  t,  this  corresponds  to  minimizing  and,  in  the  Gaussian  case 
p=2,  to  the  method  of  nonlinear  least  squares.  To  determine  the  best-fitting  x  and  t,  GSEL 
employs  a  recursive  grid-search  technique. 


Linear  Conjugate  Gradient  Inversion 

Once  the  database  of  events  has  been  relocated  using  the  appropriate  travel  time  tables, 
and  the  ray  path  sensitivities  have  been  transformed  back  onto  a  geographic  grid,  we 
perform  a  linear  conjugate  gradients  inversion  for  an  optimized  update  to  the  velocity 
model.  We  use  a  version  of  the  LSQR  algorithm  (Nolet,  1983;  Paige  and  Saunders, 

1982)  to  produce  this  update.  The  LSQR  algorithm  is  a  linear  conjugate  gradient  method 
used  to  iteratively  solve  large  systems  of  sparse,  linear  equations.  The  output  of  the 
algorithm  is  a  vector  of  changes  to  the  input  velocity  model.  The  model  update  produced 
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by  LSQR  is  constrained  in  several  ways.  First,  we  fix  a  small  buffer  region  on  the  model 
region  perimeter  to  the  values  of  the  initial  model.  This  is  to  prevent  large  velocity 
variations  from  occurring  in  areas  of  the  model  that  are  not  well  covered  by  data  and  to 
allow  us  to  seamlessly  integrate  our  final  models  into  other  global  models.  Second,  we 
apply  a  smoothing  operator  L  (from  Equation  2)  to  the  model  using  a  second  differencing 
operator,  which  is  equivalent  to  ensuring  the  curvature  of  the  model  is  smooth  (Twomey, 
1977).  Finally,  we  apply  a  scalar  damping  parameter  to  the  model  to  balance  the 
sharpness  or  noisiness  with  the  horizontal  spread  or  smoothness  of  the  recovered  velocity 
contrasts. 

After  the  linear  conjugate  gradient  method  has  converged  to  an  optimized  update  to  the 
model,  we  use  the  model  change  vector  as  a  search  direction  in  the  next  iteration  of  the 
NLCG  inversion. 


APPLICATION  TO  DATA  FROM  THE  INDIA-PAKISTAN  REGION 

In  this  section  we  detail  the  application  of  our  3-D  joint  tomography  and  location 
algorithm  to  data  from  the  India-Pakistan  region.  Figure  4  illustrates  the  setting  of  the 
WINPAK3D  velocity  model.  The  location  of  seismic  station  NIL  at  Nilore,  Pakistan 
(future  site  of  IMS  station  PRPK)  is  shown  as  the  Mack  triangle,  and  locations  of  the 
India  and  Pakistan  nuclear  test  rites  are  shown  as  blade  stars.  The  area  outlined  with  a 
black  square  is  the  contracted  study  region;  the  red  rectangle  outlines  the  region  where 
we  collected  data  to  ensure  sufficient  data  density  for  tomography  inside  the  black 
square.  Our  study  region,  centered  on  Pakistan  and  extending  into  eastern  Iran,  the 
southern  states  of  the  Former  Soviet  Union,  and  western  India,  has  a  complex  tectonic 
history  that  is  exhibited  in  the  diverse  geometries  of  crustal  structures  across  the  area. 
Mountain  ranges  extend  from  the  Kirthar  Range  in  southern  Pakistan  across  the  Sulaiman 
Range  of  central  Pakistan  and  continue  eastward  across  the  Salt  Range  in  northern 
Pakistan  and  into  the  western  Himalayas  in  India.  These  ranges  represent  a  diffuse  zone 
of  deformation  that  is  the  result  of  the  oblique  continent-continent  collision  between  India 
and  Eurasia  (Bernard  et  al.,  2000).  Further  complicating  this  deformation  zone,  a  string 
of  continental  blocks,  microcontinents,  and  island  arcs  have  been  incorporated  along  the 
Eurasian  plate  boundary  (e.g.,  Powell,  1979;  Klootwijk  etal. ,  1981;  Cobbold  and  Davy, 

1 988;  Haq  and  Davis,  1 997).  This  deformation  zone  is  bounded  to  the  west  by  the 
Chaman  fault,  a  large  sinistral  strike-slip  fault. 

The  oceanic  crust  of  the  Arabian  Sea  plate  contrasts  the  high  mountain  ranges  and  thick 
continental  crust  that  characterize  many  areas  within  the  India-Pakistan  region.  As  a 
result,  crustal  thicknesses  in  this  region  vary  widely  from  15  km  near  the  Arabian  Sea 
(Byrne  et  cd.,  1992)  to  more  than  70  km  under  the  Hindu  Kush  region  (e.g.,  Brandon  and 
Romanowicz,  1986;  Fan  et  al.,  1994).  In  addition,  there  are  several  basin  regions  where 
thick  sediments  overlie  much  thinner  continental  crust  than  is  typical  of  other  parts  of  the 
region.  The  heterogeneous  mix  of  crustal  type,  the  variability  of  crustal  thickness,  and  the 
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complex  tectonic  deformation  make  it  difficult  to  accurately  predict  travel  times  in  this 
region  without  a  realistic  regional  velocity  model. 
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Figure  4:  Physiographic  imp  of  the  region  of  interest.  The  black  rectangle  outlines  the 
area  where  the  tomographically-derived  velocity  model  is  of  the  highest  resolution;  the 
red  rectangle  outlines  the  area  where  we  collected  data  to  ensure  sufficient  ray  coverage 
in  the  tomography. 

Initial  Model 


We  compiled  a  detailed  initial  3-D  velocity  model  for  the  India-Pakistan  region  by 
synthesizing  pertinent  data  from  approximately  sixty  published  references  on  the  velocity 
structure,  geology  and  tectonics  throughout  the  region  (Figure  4).  The  references  we 
utilized  included  data  such  as  seismic  reflection  and  refraction  surveys  (i.e.  DSS  profiles, 
Pn  tomography,  Pnl  waveform  inversion),  interpretations  of  gravity  data,  surface  wave 
studies,  and  receiver  function  analyses.  Because  these  data  sources  vary  in  spatial 
coverage,  resolution,  and  the  number  of  constraints,  the  model  varies  in  a  similar  manner. 

'  The  velocity  model  is  defined  on  a  grid  of  one-degree  by  one-degree  blocks  and  5  km 
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depth  intervals  from  0  to  75  km.  We  appended  the  IASP91  model  (Kennett  and  Engdahl, 
1991)  to  the  base  of  the  preliminary  velocity  model,  beginning  at  80  km  depth  and 
extending  to  700  km  depth,  to  accommodate  ray  paths  that  travel  into  the  upper  mantle. 
Figure  5  shows  depth  slices  of  our  initial  velocity  model  at  0, 20,  40  and  60  kms 
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Figure  5:  Depth  slices  at  0  bn,  20  km,  40  km,  ami  60  km  through  the  initial  model  for 
the  India-Pabstan  model  There  is  considerable  variability  in  the  crust  and  upper  mantle 
region. 
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depth  showing  the  variability  of  velocity  model  resolution,  as  well  as  the  range  of  crustal 
thickness  across  our  study  region.  Some  preliminary  validation  was  performed  on  the 
initial  model  to  verify  its  suitability  and  potential  to  improve  event  locations.  See  Johnson 
and  Vincent  (2001)  for  details  on  the  development  and  validation  of  this  initial 
WENPAK3D  model. 

Earthquake  Dataset 

Our  data  set  is  made  up  of  events  selected  from  the  Engdahl  et  al.  (1998)  database  (EHB) 
of  well-located  earthquakes  and  explosions.  The  EHB  data  set  is  considered  to  be  GT15 
or  better  in  continental  areas,  according  to  the  IASPEI  Working  Group  on  Reference 
Events  Oittp://lemond.colorado.edu/~copgte/Y  We  chose  our  particular  subset  of  events 
from  the  EHB  database  based  on  their  spatial  distribution  across  the  region,  both  in 
latitude/longitude  and  depth.  Only  events  with  6  or  more  regional  arrivals  were  selected 
from  the  database  to  insure  sufficient  data  during  the  hypocenter  location  phase  of  the 
inversion.  In  addition,  we  selected  only  those  arrivals  with  residuals  smaller  than  7 
seconds  (i.e.  greater  than  -7  and  less  than  +7  seconds),  in  an  effort  to  filter  out  data  with 
potentially  bad  arrival  time  picks  or  incorrect  phase  assignments.  The  value  of  7  seconds 
was  chosen  to  optimize  maximum  retention  of  data,  while  still  removing  the  bulk  of  the 
poor  arrival  readings.  We  did  not  lower  the  residual  cut-off  too  much  because  residuals 
of  5  seconds  or  greater  are  possible  in  parts  of  the  model  such  as  the  Hindu  Kush  region. 
Figure  6  shows  the  histogram  of  the  residuals  in  our  data  set;  these  consist  of  the 
observed  arrival  times  minus  1ASP91  predicted  arrival  times.  By  choosing  a  cut-off  of  7 
seconds,  we  retain  96.5%  of  the  data  while  removing  those  data  with  abnormally  high 
residuals.  The  events  in  the  data  set  were  recorded  at  70  stations  within  our  model 
region,  and  consisted  of 6,626  arrivals  that  were  used  in  the  joint  inversion.  Figure  7 
depicts  the  ray  coverage  of  our  data  set,  color  coded  by  epicentral  distance  (A).  We  show 
the  data  set  divided  into  bins  of  5°  to  demonstrate  that  the  rays  to  be  used  in  the 
tomographic  inversion  are  of  the  appropriate  epicentral  distance  to  adequately  sample  the 
Pn  velocity. 

Histogram  of  Residuals  In  EHB  Dataset  for  Wln0ak3D  Model 


Figure  6:  Histogram  of  residuals  in  the  tomography  data  set.  We  retained  data  that  had 
residuals  of  7  seconds  or  less  to  reduce  the  number  of  incorrect  arrival  readings. 
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Figure  7:  Ray  coverage  (straight  line  approximations)  of  the  earthquake  data  set  used  to 
produce  a  tomographic  update  to  the  Pn  velocity  over  the  model  region.  The  ray  colors 
correspond  to  the  distance  between  station  and  event,  and  demonstrate  the  usefulness  of 
the  data  set  to  sample  Pn  adequately. 


Resolution  Test 

Prior  to  inverting  for  the  India-Pakistan  Pn  velocity  map,  we  conducted  a  resolution  test 
of  the  data  set  to  estimate  the  ability  of  the  algorithm  to  recover  a  checkerboard  pattern 
using  our  station/event  geometry.  We  overlaid  a  checkerboard  model  (shown  on  the  left 
in  Figure  8)  across  the  region  and  produced  3-D  synthetic  travel  times  for  all  the  station- 
event  pairs  in  our  earthquake  data  set.  Each  square  in  the  checkerboard  model  has  a  size 
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of  3°  in  latitude  and  3°  in  longitude.  Alternating  checkers  have  constant  velocities  of  7.9 
km/s  and  8.3  km/s,  with  a  transition  border  set  to  the  midpoint  between  high  and  low 
velocities  between  all  checkers.  Using  the  synthetic  travel  times  and  a  constant  velocity 
initial  model  (the  IASP91  Pn  velocity  of  8.04  km/s),  we  performed  one  iteration  of  our 
nonlinear  conjugate  gradient  scheme  to  retrieve  an  estimate  of  the  resolving  capability  of 
our  data  set.  The  damping  parameter  was  chosen  to  reduce  the  rms  while  keeping  the 
noise  (one-node  variations)  low.  The  damping  parameter  chosen  also  preserves  the 
amplitude  of  the  velocity  variations  of  the  synthetic  model.  The  pattern  recovered  from 
the  inversion  is  on  the  right  in  Figure  8,  with  our  specific  region  of  interest  outlined  in 
black  (between  25-40°N,  60-75°E).  Resolution  is  excellent  in  the  outlined  region,  while 
less  well-resolved  regions  to  the  south  and  northwest  reflect  the  reduced  data  coverage. 
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Figure  8:  Results  from  performing  a  checkboard  resolution  test  of  our  earthquake  data 
set  On  the  left  is  the  model  we  used  to  predict  synthetic  travel  times  from  events  to 
stations  in  our  data  set  On  the  right  is  the  checkerboard  pattern  resolved  by  one 
iteration  of  the  nonlinear  inversion  method.  Inside  the  region  of  interest  (outlined  in 
black),  the  resolution  of  the  individual  checkers  is  excellent. 
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Pn  Inversion  Results 


After  the  checkerboard  resolution  test,  we  performed  several  iterations  of  our  nonlinear 
inversion  technique  to  invert  for  an  update  to  the  Pn  velocity  map  extracted  from  our 
initial  model.  The  event  hypocenters  were  constrained  to  be  within  1 5  km  of  the  EHB 
locations,  to  correspond  with  their  GT1S  designation.  The  damping  parameter  for  the 
inversions  was  again  chosen  to  reduce  the  rms  while  keeping  the  noisiness  of  the 
recovered  velocity  change  low.  Figure  9  shows  the  initial  Pn  map  on  the  left  and  the 
model  retrieved  by  the  inversion  procedure  on  the  right.  The  general  distribution  of 
lower  and  higher  velocities  in  the  final  model  is  similar  to  the  starting  model,  but  contains 
more  detail.  The  RMS  fit  to  the  data  of  the  final  model  is  1.81,  which  is  a  significant, 
improvement  over  both  the  IASP91  value  of  2.42  and  the  starting  model  value  of  2.01. 
Table  1  lists  the  RMS  values  from  the  IASP91  model  (for  comparison)  and  each  iteration 
of  the  nonlinear  algorithm. 
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Figure  9:  Initial  (left)  and  final  (right)  Pn  velocity  maps.  The  final  Pn  velocity  model  was 
the  result  of  three  iterations  of  our  nonlinear  tomography  scheme. 
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Iteration  # 

RMS 

Reduction  in  RMS 

IASP91 

2.42 

N/A 

WO 

2.01 

0.41 

W1 

1.92 

0.09 

W2 

1.85 

0.07 

W3 

1.82 

0.03 

Table  1 :  The  RMS  results  from  several  iterations  of  the  joint  nonlinear  inversion  method 
The  first  row  entry  is  provided  for  comparison  with  the  1-D IASP91  global  model. 


In  the  region  contained  within  the  black  rectangle,  we  feel  that  the  data  coverage  is 
excellent  and  the  Pn  velocity  changes  are  reliable.  However,  one  area  of  the  model  that 
cannot  be  considered  reliable  is  to  the  southeast  on  the  continent  of  India.  This  region 
seems  well  resolved  on  the  checkerboard  model  (Figure  8),  but  is  overestimated  in  the 
updated  velocity  map  (Priestley  et  al,  2001).  High  velocities  in  this  region  can  be 
explained  by  both  poor  data  coverage  (see  the  ray  coverage  picture)  and  by  the  fact  that 
our  real  data  have  errors,  while  synthetic  data  do  not.  We  are  currently  acquiring 
additional  seismic  data  in  the  region  and  plan  to  reevaluate  the  results  in  the  future. 

It  is  somewhat  surprising  that  the  final  model  was  reached  in  only  three  iterations  of  the 
nonlinear  technique,  given  the  scale  and  underdetermined  nature  of  the  tomography 
problem.  One  possible  explanation  for  this  is  the  high  quality  of  the  initial  model,  which 
has  been  shown  to  improve  seismic  event  locations  in  the  region.  Also,  it  is  probable  that 
inverting  strictly  for  the  Pn  velocity  without  allowing  for  changes  in  the  crustal  velocities 
and  Moho  interface  depth  is  hampering  further  decreases  in  the  residuals.  Finally,  it  is 
also  possible  that  unknown  picking  errors  exist  in  the  data  that  preclude  further 
improvement  in  the  velocity  model.  Without  having  the  waveform  data  available  to 
repick  the  phases,  we  cannot  estimate  the  effect  these  errors  have  on  the  inversion 
process  or  final  model. 

After  the  third  iteration  of  the  nonlinear  inversion  method,  we  reinserted  the  final  update 
for  the  Pn  velocity  map  into  our  3-D  velocity  model  for  the  Pakistan  region.  We  then 
performed  some  preliminary  validation  on  the  model  to  verify  that  it  improves  regional 
seismic  event  location. 


Preliminary  Model  Validation 

There  is  very  little  ground  truth  data  currently  available  in  our  region  that  can  be  used  to 
validate  our  3-D  model.  However,  on  14  February  1977,  a  magnitude  5.2  earthquake 
occurred  in  the  region  near  Nilore,  Pakistan  that  was  well  recorded  by  the  Tarbela  and 
Chasma  local  networks.  Based  on  data  from  these  networks,  Seeber  and  Armbruster 
(SA)  (1979)  found  a  hypocenter  for  the  event  and  conducted  a  detailed  study  of  the 
aftershock  sequence.  Using  teleseismic  as  well  as  regional  data,  both  the  ISC  and  the 
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USGS  (NEIC)  located  the  event  within  about  5  km  of  the  epicenter  reported  by  SA,  who 
used  only  the  local  network  data.  The  ISC  calculated  the  depth  of  the  event  to  be  27  km, 
and  the  USGS  fixed  the  depth  at  33  km.  However,  the  hypocenters  of  the  main  shock 
(14.5  km)  and  the  50  accurately  located  aftershocks  determined  by  SA  indicate  a  rupture 
surface  between  12  and  18  km  depth. 

Because  of  the  regional  station  coverage  and  the  further  constraint  on  the  hypocentral 
region  derived  from  the  aftershock  study,  the  14  February  1977  event  (denoted  the 
VDAY  event)  has  one  of  the  best  constrained  locations  in  the  region.  Therefore,  we 
began  preliminary  testing  of  the  location  capability  of  our  updated  velocity  model  using  a 
subset  of  regional  stations  from  this  event.  Figure  10  shows  the  regional  stations  that 
were  used  to  locate  the  VDAY  event;  the  maximum  azimuthal  gap  is  140°.  We  used  3-D 
travel  time  tables  from  both  our  final  WINPAK3D  model  and  the  1-D IASP91  model  in 
GSEL  to  locate  the  VDAY  event. 


Figure  10:  Eleven  regional  stations  used  to  locate  the  VDAY  event  using  both  the  IASP91 
model  and  the  initial  and final  WINPAK3D  models. 


21 


Figure  1 1  shows  the  hypocenters  from  the  VDAY  event,  calculated  using  data  from  these 
1 1  regional  stations,  for  the  updated  WINPAK3D  model  (blue  x),  initial  WINPAK3D 
model  (red  x),  and  the  1-D IASP91  model  (green  x).  We  compare  these  solutions  with 
the  SA  location  for  this  event  (black  star).  The  surface  and  depth  projections  of  the 
respective  three-dimensional  confidence  regions  determined  by  Monte  Carlo  simulation 
(Rodi  and  Toksdz,  2000)  show  the  95%  confidence  levels  for  each  model's  hypocenter. 
The  epicenter  mislocation  of  the  updated  WINPAK3D  model  from  the  SA  event  location 
is  3.5  km,  while  the  initial  model's  epicenter  mislocation  is  9.3  km  and  the  IASP91 
epicenter  mislocation  is  31.7  km.  In  addition,  the  95%  confidence  region  for  both  the 
initial  and  updated  WINPAK3D  model  epicenters  encompass  the  SA  epicenter  while  the 
95%  confidence  region  for  the  IASP91  epicenter  does  not. 

Also  shown  in  Figurel  1  are  the  hypocenters  determined  by  both  the  ISC  and  the  USGS 
(NEIC)  (open  circles)  using  leleseismic  data  as  wed!  as  regional  data.  Note  that  by  using 
regional  data  alone,  foe  3-D  model  is  able  to  estimate  foe  hypocenter  of  this  event  as  well 
as  foe  teleseismically-derived  estimates.  These  results  are  encouraging  for  the  potential 
success  of  using  3-D  regional  velocity  models  to  find  accurate  locations  of  small  events 
that  are  not  recorded  teleseismicaUy. 


Figure  11:  Hypocenters from  the  VDA  Y  event,  calculated  using  data  from  11  regional 
stations,  for  the  updated  WINPAK3D  model  (blue  x),  initial  WINPAK3D  model  (red  x), 
cmd  the  1-D  IASP91  model  (green  x).  We  compare  these  solutions  with  the  SA  location 
for  this  event  (black  star).  Also  shown  are  the  hypocenters  determined  by  both  the  ISC 
and  the  USGS  (NEIC)  (open  circles)  using  teleseismic  data  as  well  as  regional  data 


Figure  12  shows  a  source-specific  station  correction  (SSSC)  derived  from  the 
WINPAK3D  regional  velocity  model  for  foe  India-Pakistan  region  for  a  source  at  15  km 
depth.  This  correction  surface  was  produced  by  calculating  travel  times  from  foe  surface 
at  foe  Nilore  (NIL)  station  to  all  other  points  in  the  3-D  model,  which  was  discretized  on 
a  5  km  grid.  We  then  subtracted  foe  travel  times  through  the  IASP91  model,  which  we 
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discretized  onto  an  equivalent  5  km  grid  and  derived  using  the  identical  PL  travel  time 
prediction  method.  Figure  12  represents  the  anticipated  station  correction  with  respect  to 
the  LA.SP91  model  for  an  event  occurring  at  1 5  km  depth  beneath  station  NIL. 
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Figure  12:  Model-based  source  specific  station  correction  (SSSC)  at  station  NIL  for 
events  occurring  at  15  km  depth. 


Since  the  VD  AY  event  occurred  near  the  Nilore  station,  it  provides  a  unique  reciprocity 
test  of  our  3-D  velocity  model.  Stations  that  recorded  the  event  in  other  parts  of  our 
model  act  as  sources  that  could  potentially  propagate  seismic  energy  to  the  station  at 
Nilore.  We  analyzed  the  residuals  from  this  event  and  compared  them  with  the  SSSC 
derived  from  the  WINPAK3D  model.  In  Figure  13  we  show  this  comparison,  which 
illustrates  the  difference  between  the  empirical  travel  time  corrections  (residuals)  and  the 
corrections  derived  from  the  WINPAK3D  model.  The  two  are  in  overall  agreement,  with 
absolute  differences  generally  from  0  to  3  seconds.  Note  that  SSSCs  can  range  up  to  8 
seconds  in  some  regions.  For  approximately  half  the  stations,  we  have  a  less  than  1 
second  difference. 
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Figure  13:  Differences  between  empirical  travel  time  corrections  (residuals)  cmd  the 
corrections  derived  from  the  WINPAK3D  model  for  all  23  regional  stations  used  to 
locate  the  VDAY event. 


CONCLUSIONS  AND  RECOMMENDATIONS 

Our  primary  goal  in  this  project  was  to  develop  a  3-D  model  of  the  crust  and  upper 
mantle  in  the  Pakistan  region  to  improve  regional  seismic  event  location.  We  developed 
our  3-D  model  by  applying  a  joint  velocity  tomography  and  event  location  procedure  to  a 
dataset  of  earthquakes  from  the  region.  The  model  was  iteratively  updated  using  a 
nonlinear,  conjugate-gradients  technique  that  adjusts  the  velocity  model  to  minimize  the 
misfit  between  calculated  and  observed  travel  times  from  multiple  stations  and  events, 
subject  to  smoothness  constraints.  We  performed  several  iterations  of  this  nonlinear 
algorithm  to  invert  for  the  Pn  velocity  as  a  function  of  latitude  and  longitude  and  then 
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imposed  a  fixed  rule  for  extrapolating  this  velocity  into  the  upper  mantle.  Our  results 
show  that  the  new  model  fits  the  data  better  compared  to  both  the  initial  model  and  the 
global  1-D IASP91  model.  The  root  mean  square  (RMS)  error  for  the  updated  3-D  model 
is  1.81,  compared  to  2.01  for  the  initial  3-D  model  and  2.42  for  the  IASP91  model. 
Preliminary  validation  of  the  model  indicates  that  it  does  improve  regional  seismic  event 
location. 

This  algorithm  can  be  effectively  transported  for  use  in  other  regions  of  the  world  and  has 
already  been  applied  to  the  area  surrounding  the  IMS  station  BRVK  at  Borovoye,  Russia 
(Murphy  et  al.,  2001).  It  is  currently  one  of  the  most  powerful  methods  available  to 
produce  tomographic  3-D  regional  velocity  models.  These  models  can  then  be  used  to 
generate  source-specific  station  corrections  (SSSC’s)  for  use  in  accurate  regional  and 
teleseismic  event  locations.  We  plan  further  extensions  to  the  method,  including  allowing 
for  more  flexible  changes  in  the  upper  crust  and  mantle  portions  of  the  model  and 
improving  the  ray  tracing  algorithm. 
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